Project_V1

suppressPackageStartupMessages(library(readxl))
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(plotly))
suppressPackageStartupMessages(library(fmsb))

data = read_excel("Project_1_Data.xlsx", sheet = "pooled123")

filteredData = select(data, Date, PID, BSSQ_1:BSSQ_16, ASSQ_1:ASSQ_16, age, VRexperience, ssq_full, gender, sm)
#we only want to look at those that did not undergo social modelling conditions
filteredData = subset(filteredData, sm == "NO_SM")

#reclasss VR experience as factor (was chr)
filteredData$VRexperience = as.factor(filteredData$VRexperience)

#we want to filter this data even further and split it into age groups
#once in age groups, calculate the mean change for each of the age groups for each symptom
filteredData = mutate(filteredData, age_group = case_when(
  age > 45 ~ "above 45",
  age >= 38 ~ "38 to 45",
  age >= 30 ~ "30 to 37",
  age >= 22 ~ "22 to 29",
  age >= 16 ~ "16 to 21"
))

symptoms = c("discomfort",
                "fatigue", 
                "headache",
                "eyestrain",
                "difficulty_focusing",
                "salivation",
                "sweating",
                "nausea",
                "difficulty_concentrating",
                "fullness_of_head",
                "blurred_vision",
                "dizziness_o",
                "dizziness_c",
                "vertigo",
                "stomach_awareness",
                "burping")

for (i in 1:length(symptoms)) {
  filteredData[[paste0("d_", symptoms[[i]])]] <- filteredData[[paste0("ASSQ_", i)]] - filteredData[[paste0("BSSQ_", i)]]
} # Borna's friend made fun of us for doing this manually so we got chatGPT to fix it

yes_means_vector <- numeric(length(symptoms))
no_means_vector = numeric(length(symptoms))

# Loop through symptoms and calculate the mean for each one
for (i in 1:length(symptoms)) {
  # Calculate the mean for the current symptom column
  yes_means_vector[i] <- mean(filter(filteredData, VRexperience == "Yes")[[paste0("d_", symptoms[i])]])
  no_means_vector[i] <- mean(filter(filteredData, VRexperience == "No")[[paste0("d_", symptoms[i])]])
}

# Create a data frame from the vector
means_per_symptom <- data.frame(symptom = symptoms, yes_mean_value = yes_means_vector, no_means_value = no_means_vector)

Summary of Findings

(executive summary goes here)

Initial Data Analysis (IDA)

Source

Our data was sourced from Cosette Saunder’s paper “Socially Acquired Nocebo Effects Generalize but Are Not Attenuated by Choice”.

Structure

The paper included 336 participants, each with 51 variables recorded. Our project focused the following variables: We specifically looked at participants that did not undego any social modelling; 69 participants, each with 51 variables recorded. Our project focused the following variables:

• Baseline SSQ (BSSQ), Active SSQ (ASSQ), and SSQ_Fullof 16 symptoms (quantitative, discrete): self-reported symptom severity before and after undergoing VR respectively, on a scale of 1 to 10.

• Participants’ VRexperience (qualitative, nominal); sorted by experience to see how symptoms differed.

• age of participants (quantitative, discrete); they were then sorted into age groups.

• The change (\(\Delta\)) between BSSQ and ASSQ was calculated for each participant for each symptom, and we took the average \(\Delta\) for each symptom for each group (VRexperience and age_group) (quantitative, discrete).

Age Groups

data_age_groups = select(filteredData, age_group)
groups_counted = data_age_groups %>% count(age_group)

age_pie = plot_ly(groups_counted, labels = ~age_group, values = ~n,
                 type = 'pie', direction = 'clockwise', sort = FALSE, rotation = 30)
age_pie <- age_pie %>% layout(title = 'Distribution of ages',
                            showlegend = TRUE)
age_pie

VR Experience

data_experience = select(filteredData, VRexperience)
exp_counted = data_experience %>% count(VRexperience)

vr_pie = plot_ly(exp_counted, labels = ~VRexperience, values = ~n,
                 type = 'pie')
vr_pie <- vr_pie %>% layout(title = 'Distribution of VR experience',
                            showlegend = TRUE)

vr_pie

Limitations

• Self-reporting bias for both BSSQ and ASSQ which makes the values prone to being over or underestimates by each participant.

• VRexperience is a binary qualitative classification, it is not descriptive of the nature, amount or frequency, and that those with more than 10 VR experiences were excluded.

Research Question

Research Question

What is the effect of having VR experience on the symptoms experienced by people?

What is the effect of past VR experience on the symptoms experienced by people?

plt = ggplot(filteredData, aes(x = VRexperience, y = ssq_full, fill = VRexperience)) +
  geom_boxplot() +
  theme_minimal() +
  labs(x = "VR Experience", y = "Full SSQ")
ggplotly(plt)

Those with and without VR experience reported varied severity of symptoms. Interestingly, the VR experienced group saw a greater median ssq_full, at 7.0 compared to -0.5 for those without experience, indicating those with experience reported more severe symptoms afterwards. VR experienced participants also showed a greater IQR than those without, at 16.25 and 8.5 respectively, pointing to experienced users reporting less consistent symptom severity.

library(plotly)

fig <- plot_ly(
  type  = 'scatterpolar',
  fill = 'toself',
  title = "Average change in BSSQ and ASSQ depending on VR experience for each symptom"
)

fig <- fig %>% 
  add_trace(
    r = means_per_symptom$yes_mean_value,
    theta = c('Discomfort', 'Fatigue', 'Headache', 'Eyestrain', 'Sweating', 'Nausea', 'Dizziness (closed)', 'Dizziness (open)', 'Difficulty Focusing', 'Difficulty Concentrating', 'Salivation', 'Blurry Vision', 'Vertigo', 'Stomach Awareness', 'Fullness of head','Burping'),
    name = 'With VR Experience'
  )

fig <- fig %>% 
  add_trace(
    r = means_per_symptom$no_means_value,
    theta = c('Discomfort', 'Fatigue', 'Headache', 'Eyestrain', 'Sweating', 'Nausea','Dizziness (closed)', 'Dizziness (open)', 'Difficulty Focusing', 'Difficulty Concentrating', 'Salivation', 'Blurry Vision',  'Vertigo', 'Stomach Awareness', 'Fullness of head','Burping'),
    name = 'No VR Experience'
  )

fig <- fig %>%
  layout(
    polar = list(
      radialaxis = list(
        visible = T,
        range = c(-0.5,2)
      )
    ),
    showlegend = T
  )

fig

When took the average change between the BSSQ and ASSQ, and visualized it in the spiderchart. We can see that both groups experience similar \(\Delta\) (change) in most symptoms. Interestingly, the VR experienced group saw The spider-chart above reinforces what we see in the box-plot, with VR experienced users all indicating a higher average \(\Delta\) for almost all of the 16 symptoms.

This observation is in contrast to These observations are in contrast to [insert article(s) here]

We look further into those symptoms with the greatest differences in mean \(\Delta\) : Dizziness (Eyes closed) , Vertigo and Stomach Awareness.

ggplot(filteredData) +
  geom_boxplot(aes(x = "Discomfort", y = d_discomfort, fill = VRexperience)) +
  geom_boxplot(aes(x = "Fatigue", y = d_fatigue, fill = VRexperience)) +
  geom_boxplot(aes(x = "Headache", y = d_headache, fill = VRexperience)) +
  geom_boxplot(aes(x = "Eyestrain", y = d_eyestrain, fill = VRexperience)) +
  geom_boxplot(aes(x = "Difficulty Focusing", y = d_difficulty_focusing, fill = VRexperience)) +
  geom_boxplot(aes(x = "Salivation", y = d_salivation, fill = VRexperience)) +
  geom_boxplot(aes(x = "Sweating", y = d_sweating, fill = VRexperience)) +
  geom_boxplot(aes(x = "Nausea", y = d_nausea, fill = VRexperience)) +
  geom_boxplot(aes(x = "Difficulty Concentrating", y = d_difficulty_concentrating, fill = VRexperience)) +
  geom_boxplot(aes(x = "Fullness of Head", y = d_fullness_of_head, fill = VRexperience)) +
  geom_boxplot(aes(x = "Blurred Vision", y = d_blurred_vision, fill = VRexperience)) +
  geom_boxplot(aes(x = "Dizziness (o)", y = d_dizziness_o, fill = VRexperience)) +
  geom_boxplot(aes(x = "Dizziness (c)", y = d_dizziness_c, fill = VRexperience)) +
  geom_boxplot(aes(x = "Vertigo", y = d_vertigo, fill = VRexperience)) +
  geom_boxplot(aes(x = "Stomach Awareness", y = d_stomach_awareness, fill = VRexperience)) +
  theme_minimal() +
  labs(title = "VR Experience and Symptom Change",
       x = "Symptom",
       y = "Change in Severity") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

library(tidyverse)
library(plotly)

plt = ggplot(filteredData, aes(x = as.factor(gender) , y = ssq_full, group = VRexperience, fill = VRexperience)) +
  geom_boxplot(position = "dodge") +
  labs(x = "gender",  y = "Full SSQ")
plt

ggplotly(plt) %>%  layout(boxmode = "group")

Warning: 'layout' objects don't have these attributes: 'boxmode'
Valid attributes include:
'_deprecated', 'activeshape', 'annotations', 'autosize', 'autotypenumbers', 'calendar', 'clickmode', 'coloraxis', 'colorscale', 'colorway', 'computed', 'datarevision', 'dragmode', 'editrevision', 'editType', 'font', 'geo', 'grid', 'height', 'hidesources', 'hoverdistance', 'hoverlabel', 'hovermode', 'images', 'legend', 'mapbox', 'margin', 'meta', 'metasrc', 'modebar', 'newshape', 'paper_bgcolor', 'plot_bgcolor', 'polar', 'scene', 'selectdirection', 'selectionrevision', 'separators', 'shapes', 'showlegend', 'sliders', 'smith', 'spikedistance', 'template', 'ternary', 'title', 'transition', 'uirevision', 'uniformtext', 'updatemenus', 'width', 'xaxis', 'yaxis', 'barmode', 'bargap', 'mapType'

Dizziness

library(tidyverse)
library(plotly)

plt = ggplot(filteredData, aes(x = VRexperience, y = d_dizziness_c)) +
  geom_boxplot() +
  labs(x = "VR Experience",  y = "Change in BSSQ and ASSQ-Dizziness (Eyes Closed)")

ggplotly(plt)

#fivenum(withVRexperience$d_dizziness_c)
#fivenum(noVRexperience$d_dizziness_c)
#
#mean(withVRexperience$d_dizziness_c)
#mean(noVRexperience$d_dizziness_c)
#
#sd(withVRexperience$d_dizziness_c)
#sd(noVRexperience$d_dizziness_c)


#since from the boxplots we can see that both sets have a fair number of outliers which could affect the mean, the median + IQR may be a better measure of spread

library(tidyverse)
library(plotly)
library(gganimate)
library(gifski)

plt = ggplot(filteredData, aes(y = d_dizziness_c, x = age, colour = VRexperience)) +
  geom_point() +
  labs(x = "Age",  y = "Dizziness (Eyes Closed)", col = "VR Experience") +
  transition_states(VRexperience, transition_length = 1, state_length = 1) +
  enter_fade() + exit_fade()

#animate(plt, renderer = gifski_renderer())
ggplotly(plt)

library(tidyverse)
library(plotly)

plt = ggplot(filteredData, aes(x = age_group, y = d_dizziness_c)) +
  geom_boxplot() +
  labs(x = "Age Group",  y = "Change in BSSQ and ASSQ-Dizziness(Eyes Closed)")

ggplotly(plt)

Vertigo

library(tidyverse)
library(plotly)

plt = ggplot(filteredData, aes(x = VRexperience, y = d_vertigo)) +
  geom_boxplot() +
  labs(x = "VR Experience",  y = "Vertigo")

ggplotly(plt)

#fivenum(withVRexperience$d_vertigo)
#fivenum(noVRexperience$d_vertigo)
#
#mean(withVRexperience$d_vertigo)
#mean(noVRexperience$d_vertigo)
#
#sd(withVRexperience$d_vertigo)
#sd(noVRexperience$d_vertigo)

library(tidyverse)
library(plotly)

plt = ggplot(filteredData, aes(x = age_group, y = d_vertigo)) +
  geom_boxplot() +
  labs(x = "Age Group",  y = "Change in BSSQ and ASSQ - Vertigo")

ggplotly(plt)

library(tidyverse)
library(plotly)
library(gganimate)
library(gifski)

plt = ggplot(filteredData, aes(y = d_vertigo, x = age, colour = VRexperience)) +
  geom_point() +
  labs(x = "Age",  y = "Vertigo", col = "VR Experience") +
  transition_states(VRexperience, transition_length = 1, state_length = 1) +
  enter_fade() + exit_fade()

#animate(plt, renderer = gifski_renderer())
ggplotly(plt)

Stomach Awareness

library(tidyverse)
library(plotly)

plt = ggplot(filteredData, aes(x = VRexperience, y = d_stomach_awareness)) +
  geom_boxplot() +
  labs(x = "VR Experience",  y = "Change in BSSQ and ASSQ - Stomach Awareness")

ggplotly(plt)

library(tidyverse)
library(plotly)

plt = ggplot(filteredData, aes(x = age_group, y = d_stomach_awareness)) +
  geom_boxplot() +
  labs(x = "Age Group",  y = "Change in BSSQ and ASSQ - Stomach Awareness")

ggplotly(plt)

library(tidyverse)
library(plotly)
library(gganimate)
library(gifski)

plt = ggplot(filteredData, aes(y = d_stomach_awareness, x = age, colour = VRexperience)) +
  geom_point() +
  labs(x = "Age",  y = "Stomach Awareness", col = "VR Experience") +
  transition_states(VRexperience, transition_length = 1, state_length = 1) +
  enter_fade() + exit_fade()

#animate(plt, renderer = gifski_renderer())
ggplotly(plt)

library(tidyverse)

p = ggplot(filteredData, aes(x = age, y = ssq_full, colour = VRexperience)) +
  geom_point()

p

plt = ggplot(filteredData, aes(x = VRexperience, y = ssq_full)) +
  geom_boxplot()
plt = ggplot(filteredData, aes(x = VRexperience, y = ssq_full, fill = VRexperience)) +
  geom_boxplot() +
  theme_minimal()

ggplotly(plt)

plt2 = ggplot(filteredData, aes(x = age_group, y = ssq_full)) +
  geom_boxplot()

ggplotly(plt2)

library(tidyverse)

yex = filter(filteredData, VRexperience == "Yes")
nox = filter(filteredData, VRexperience == "No")

q1 = quantile(yex$age, probs = c(0.25), na.rm = TRUE) |> as.numeric()
q2 = quantile(nox$age, probs = c(0.25), na.rm = TRUE) |> as.numeric()
# https://forum.posit.co/t/make-the-r-function-quantile-return-a-numeric-value/178084
q3 = 2*q2-q1

ggplot(filteredData, aes(x = VRexperience, y = age)) + geom_boxplot() + geom_hline(yintercept = q1) + geom_hline(yintercept = q2) + geom_hline(yintercept = q3)

filteredData = mutate(filteredData, age_class = case_when(age > q3 ~ "older",
                                                      age > q2 ~ "old",
                                                      age > q1 ~ "young",
                                                      age > 0 ~ "younger"))

ggplot(filter(filteredData, age_class == "old" | age_class == "young"), aes(x = age_class, y = ssq_full)) + geom_boxplot()

filteredData["age_group"] = cut(filteredData$age, c(16, 22, 30, 38, 46, Inf), c("16-21", "22-29", "30-37", "38-45", "45+"), include.lowest = TRUE)

ggplot(filteredData, aes(x = VRexperience, y = ssq_full, fill = VRexperience)) + stat_summary(fun = "mean", geom = "bar", position = "dodge")

ggplot(filteredData) +
  geom_boxplot(aes(x = "Discomfort", y = d_discomfort, fill = VRexperience)) +
  geom_boxplot(aes(x = "Fatigue", y = d_fatigue, fill = VRexperience)) +
  geom_boxplot(aes(x = "Headache", y = d_headache, fill = VRexperience)) +
  geom_boxplot(aes(x = "Eyestrain", y = d_eyestrain, fill = VRexperience)) +
  geom_boxplot(aes(x = "Difficulty Focusing", y = d_difficulty_focusing, fill = VRexperience)) +
  geom_boxplot(aes(x = "Salivation", y = d_salivation, fill = VRexperience)) +
  geom_boxplot(aes(x = "Sweating", y = d_sweating, fill = VRexperience)) +
  geom_boxplot(aes(x = "Nausea", y = d_nausea, fill = VRexperience)) +
  geom_boxplot(aes(x = "Difficulty Concentrating", y = d_difficulty_concentrating, fill = VRexperience)) +
  geom_boxplot(aes(x = "Fullness of Head", y = d_fullness_of_head, fill = VRexperience)) +
  geom_boxplot(aes(x = "Blurred Vision", y = d_blurred_vision, fill = VRexperience)) +
  geom_boxplot(aes(x = "Dizziness (o)", y = d_dizziness_o, fill = VRexperience)) +
  geom_boxplot(aes(x = "Dizziness (c)", y = d_dizziness_c, fill = VRexperience)) +
  geom_boxplot(aes(x = "Vertigo", y = d_vertigo, fill = VRexperience)) +
  geom_boxplot(aes(x = "Stomach Awareness", y = d_stomach_awareness, fill = VRexperience)) +
  theme_minimal() +
  labs(title = "VR Experience and Symptom Change",
       x = "Symptom",
       y = "Change in Severity") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

filteredData = mutate(filteredData, age_group = case_when(age >= 38 ~ "38 to 45",
                                                          age >= 37 ~ "30 to 37",
                                                          age >= 22 ~ "22 to 29",
                                                          age >= 16 ~ "16 to 21"))

ggplot(filteredData, aes(x = age_group, y = ssq_full, fill = VRexperience)) + geom_boxplot()

Professional Standard of Report

Acknowledgements

References